Testing Environment for Eurus

In [1]:

    

import numpy as np
from eurus import Eurus


    

In [2]:

    

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib

%matplotlib inline

from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 150 # Change this to adjust figure size


    

In [3]:

    

class Source(object):
    
    def __init__(self, sc):
                
        self._x, self._z = np.mgrid[
            0:sc['dx']*sc['nx']:sc['dx'],
            0:sc['dz']*sc['nz']:sc['dz']
        ]
    
    def __call__(self, x, z):
        
        dist = np.sqrt((self._x - x)**2 + (self._z - z)**2)
        srcterm = 1.*(dist == dist.min())
        nz, nx = self._x.shape
        
        return np.hstack([srcterm.T.ravel() / srcterm.sum(), np.zeros(nx*nz, dtype=np.complex128)])


    

In [4]:

    

# Geometry parameters
dx          = 1.
dz          = 1.
nx          = 100
nz          = 200

# Bulk parameters
velocity    = 2000.     * np.ones((nz,nx))
density     = 1.        * np.ones((nz,nx))

# Anisotropy parameters
theta       = 0.        * np.ones((nz,nx))
epsilon     = 0.2       * np.ones((nz,nx))
delta       = 0.2       * np.ones((nz,nx))

# Other parameters
freq        = 200.
freeSurf    = [False, False, False, False]
nPML        = 10


# Pack values into systemConfig dictionary
systemConfig = {
    'nx':       nx,
    'nz':       nz,
    'dx':       dx,
    'dz':       dz,

    'c':        velocity,
    'rho':      density,
    
    'theta':    theta,
    'eps':      epsilon,
    'delta':    delta,

    'freq':     freq,
    'freeSurf': freeSurf,
    'nPML':     nPML,
    'cPML':     1e3,
}


    

In [5]:

    

Ainv = Eurus(systemConfig)
src = Source(systemConfig)


    

In [6]:

    

plt.spy(Ainv.A, markersize=0.1, marker=',')


    

    
Out[6]:





<matplotlib.lines.Line2D at 0x10745f450>






    

In [7]:

    

q = src(50,100)
u = Ainv*q


    

In [8]:

    

clip = 1e-1

plotopts = {
    'extent':   [0, nx*dx, nz*dz, 0],
    'cmap':     cm.bwr,
    'vmin':     -clip,
    'vmax':     clip,
}

fig = plt.figure()

ax = fig.add_subplot(1,1,1, aspect=1)
plt.imshow(u[:nx*nz].reshape((nz,nx)).real, **plotopts)
plt.title('Real Wavefield')
plt.xlabel('X')
plt.ylabel('Z')


    

    
Out[8]:





<matplotlib.text.Text at 0x107440b90>






    

In [ ]:

    




    

In [ ]: